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Abstract 

In addition to their self-renewal capabilities, hematopoietic stem cells guarantee 
the continuous supply of fully differentiated, functional cells of various types in 
the peripheral blood. The process which controls differentiation into the different 
lineages of the hematopoietic system (erythroid, myeloid, lymphoid) is referred to 
as lineage specification. It requires a potentially multi-step decision sequence which 
determines the fate of the cells and their successors. It is generally accepted that 
lineage specification is regulated by a complex system of interacting transcription 
factors. However, the underlying principles controlling this regulation are currently 
unknown. 

Here, we propose a simple quantitative model describing the interaction of two 
transcription factors. This model is motivated by experimental observations on 
the transcription factors GATA-1 and PU.l, both known to act as key regulators 
and potential antagonists in the erythroid vs. myeloid differentiation processes of 
hematopoietic progenitor cells. We demonstrate the ability of the model to account 
for the observed switching behavior of a transition from a state of low expression of 
both factors (undifferentiated state) to the dominance of one factor (differentiated 
state). Depending on the parameter choice, the model predicts two different pos- 
sibilities to explain the experimentally suggested, stem cell characterizing priming 
state of low level co-expression. Whereas increasing transcription rates are suffi- 
cient to induce differentiation in one scenario, an additional system perturbation 
(by stochastic fluctuations or directed impulses) of transcription factor levels is 
required in the other case. 
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1 Introduction 



The hematopoietic system consists of a variety of functionally different cell 
types, including mature cells such as erythrocytes, granulocytes, platelets, or 
lymphocytes, as well as several different pre cursor cells (i.e., p r ematu re cell 
stages ) and hematopoietic stem cells (BSC) Q, ll9Q7t lOrkinl . 1200(1) . Most 
mature cell types have limited life spans ranging from a few hours to several 
months, which implies the existence of a source capable of replenishing these 
differentiated cells throughout the life span of an individual. This supply is 
realized by the population of HSC, which is maintained and even regenerated 
after injury or depletion throughout the whole life of th e organism. This se lf- 



renew al property is a ma.ior cnaracteristicaennmg ribU (^oemer an a rioeaer , 
20021 : lLordl . fll)97t IPotten and Loefflerl . Il990h . A second major characteristic of 
HSC is their ability to contribute to the production of cells of all hematopoietic 
lineages, thus ensuring the supply of functionally differentiated cells meeting 
the needs of the organism. The process controlling the development of undif- 
ferentiated stem or progenitor cells into one specific functional direction (i.e., 
one specific hematopoietic lineage) is called lineage specification. It is gener- 
ally accepted that the process of lineage spe cification is governed by the inter- 



play of many diffe r ent transcrip t ion factors (I Akashi 



2002 ICross et all 11994 lOrkinl Il995l 1200(1 



Vkashi, 
Tenenl . 



2005; Cantor and Orkin 



2003). Experimental re- 



sults suggest that a number of relevant trans cription factors are expressed 
simultaneously in HSC, although at a low level ( Akashi et al. , 20031 : Hu et al. . 
1997). Some authors refer to this state of a low level co-expression as vrim - 



ing behavior ( Akashi . 20051: Cross and Enver . 1997 : Enver and Greavesl . 1998). 
During differentiation the balanced co-expression of these potentially antag- 
onistic transcription factors is assumed to be broken at some point (or even 
multiple points). Thereafter, the system is supposed to be characterized by 
an up-regulated level of some transcription factors, specific for a particular 
lineage, while other transcription factors are down-regulated. These observa- 
tions suggest a transcription factor network, capable of switch-like behavior 
by changing from unspecific co-expression to different states of specific expres- 
sion. However, the general underlying principles of the regulatory mechanisms 
are currently unknown. Particularly, it is unclear whether the assumption of 
a dynamically balanced low level co-expression state is justified or whether 
priming should rather be interpreted as the result of an inactive transcrip- 
tion factor network overlaid by stochastic fluctuations of transcription factor 
expression. 



In this paper we propose a simple mathematical model describing different 
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interaction scenarios of two transcription factors. Biologically, this simple 
two component network model is motivated by experimental observations 
on the transcription factors GATA-1 and PU.l, known to be involved in 



the process of lineage speci f ication of HSC ( Du et al. ._ 2002t Oikawa et al 



1999; iRekhtman et all [l999; Rosmarin et al 



200.4 iTenenl . 120031 : IVoso et al 



1994). The zinc finger factor GATA-1 is reported to be required for the dif- 
ferentiation and maturation of erythroid/megakaryocytic cells, while the Ets- 
family transcription factor PU.l supports the deve l opmen t of m yeloid and 
lymphoid cells (reviewed by Cantor and Orkin . 20021: Tenen . 20031 ). For both, 
GATA-1 and PU.l, it has been demonstrated that they are able to stimu- 



late t hdr_owji_Jrai]^ip_y ; on through an auto-catalytic process ( Chen et al. 
199.4 iNishimura et all Eoool : lOkuno et all 1200.4 fTsai et all Il99ll ). Addition" 



ally, there are physical interactions between GATA-1 and PU.l which induce a 



myeloid) | 


Du et al . 20021 Nerlov et al. 


Voso et al. 


. 1994 Yamada et al.. 1998 



mu- 



tual inhibition of these two transcription factors have been suggested by ex- 
perimental observations: On one hand, GATA-1 binds to the (33/(34 region of 
PU.l (complex 1) and displaces the PU.l co-activator c-Jun f rom its bind- 
ing s ite, thereby, inhibiting the transcription initiation of PU.l ((Zhang et al. 
1999). On the other hand, the inhibition of GATA-1 transcription is medi 



ated by the binding of the N-terminal region of PU.l to the C-finger region 
of GA TA-1 (complex 2), thus blocking the binding of GATA-1 to its pro- 
moter ( Zhang et al. . 2000j ). That means, although both inhibition mechanisms 
are interfered through the formation of PU.l/GATA-1 heterodimers, the two 
complexes are structurally different. Whereas complex 1 (inhibition of PU.l 
transcription by GATA-1) is known to bind to DNA, thus occupying a PU.l 
promoter site, DNA-binding of complex 2 (inhibition of GATA-1 transcription 
by PU.l) has not been reported so far. 



The mechanisms of antagonistic interdependence together with positive auto- 
catalytic regulation provide a framework for the theoretical investigation of 
different scenarios of transcription factor interaction and their implications 
for the explanation of lineage specification control. Applying a mathematical 
model, which formalizes the described interactions, it is now possible to ana- 
lyze different combinations of transcription factor activation and inhibition on 
a qualitative and quantitative level. The proposed model relie s on principles 



suggested for the description of general genetic switches /e.g. iBecskei et al 
l200l[ ICinauin and Demongeoti bool 1200.4 iGardner et all . l2000i) . 



In this paper it is our objective to examine the following questions within the 
framework of this model structure: 



Are the experimentally described interactions of the two transcription fac- 
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tors sufficient to generate a switching behavior between a stable co-expression 
of two factors and the dominance of one of these factors? 

• What are the conditions inducing such a qualitative change in the system 
behavior? 

• Is there evidence for a functional role of the (experimentally suggested) 
priming status? 

To answer these question the following strategy is applied. Firstly, the model 
equations are derived on the basis of the described biological mechanisms of 
transcription factor interaction for GATA-1 and PU.l (Section 2). Secondly 
this model is analyzed with respect to the existence of steady state solutions 
and their dependence on the model parameters. According to our objective, 
to understand the mechanisms leading to switches between different stable 
system states, we focus our analysis particularly on the determination of bi- 
furcation conditions, considering different scenarios of transcription factor in- 
teraction (Section 3). Finally, the obtained results are discussed in relation to 
the ongoing debate about lineage specification control in the hematopoietic 
system, specifically with respect to potential explanations of the experimen- 
tally suggested low level co-expression of transcription factors (priming) in 
undifferentiated progenitors and stem cells (Section 4). 



2 Model description 

Although our analysis is motivated by experimental observations of specific 
transcription factor interactions (GATA-1 and PU.l), the model may also be 
applied in the general context of two interacting transcription factors. In the 
following, the two transcription factors are denoted by X and Y . 

2.1 General assumptions 

The general design of the model structure is based on the following assump- 
tions which are motivated by the experimental observations outlined in Section 
1: 

• Both transcription factors, X and Y, are able to act as activator molecules: 

- If bound to their own promoter region, X and Y introduce a positive 
feedback on their own transcription. This process is referred to as specific 
transcription (Fig. 1(a)). 

- X and Y are both able to induce an overall transcription which also 
effects potentially antagonistic transcription factors. Although such an 
interaction is most likely indirect, for the model we consider a mutual 
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activation of X and Y by the opposing transcription factor, which we 

refer to as unspecific transcription (Fig. 1(b)). 
We assume that transcription initiation is only achieved by the simultaneous 
binding of two X and Y molecules, respectively (i.e., binding cooperativity 
c = 2). This assumption is motivated by the result that a binding coopera- 
tivity c > 1 is a necessary condition for the existence of system bistability 
(see e .g. Becskei et al. . 2001 : Cinquin and Demongeot . 20051: Gardner et al 
2000j. 



There is a mutual inhibition of X and Y. Within this context, two possible 
mechanisms, based on the formation of two structurally different complexes 
of X and Y, are considered: 

- Joint binding of X and Y molecules to promoter sites (Fig. 1(c)). Here, 
the DNA-bound Xy-complex (Zi) acts as a transcription repressor, which 
blocks the binding sites. This represents a mode of competition for free 
binding sites. 

- Formation of another XF-complex, called Z2, which neither binds to X 
nor Y DNA binding site (Fig. 1(d)). In contrast to Zi, this represents a 
competition for free transcription factor molecules. 

Both inhibition mechanisms (including combinations of them) are consid- 
ered for X as well as for Y. 



To facilitate the analysis of the mathematical model we make the following 
simplifications: 



• Post-transcriptional regulation is neglected, i.e., the transcription of a gene 
is considered to ultimately result in the production of the corresponding 
protein (here, a transcription factor). 

• Time delays due to transcription and translation processes are neglected. 

• Simultaneous binding of X/Y monomers together with a Zi-heterodimer, 
of two Zi-heterodimers, as well as of a X and a Y monomer at the same 
promoter are excluded from the analysis. 

• Interactions of X, Y as well as the promoter regions of the coding genes 
with further transcription factors are neglected. 



Throughout the paper the following notations are used: x, y denote the molecule 
concentrations of X and Y, respectively. Z\ denotes the DNA bound XY- 
complex and Z<i the structurally different XY- complex, which is not able to 
bind to promoter DNA. D x / y denotes free DNA binding sites within the pro- 
moter region of X and Y, respectively. In contrast, binding sites occupied by 
X or Y molecules or by the XF-complex Z\ are denoted as jj^l yy / xy _ 
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Figure 1. Principles of transcription initiation and inhibition for X and Y. (a) 
Specific transcription, i.e. auto-catalysis by the transcription factor itself; (b) Un- 
specific transcription, i.e. transcription initiated by another transcription factor; 
(c-e) Suggested mechanisms of transcription inhibition for X and Y by formation 
of XY-complexes: (c) A XY-complex, called Z\, bound to the promoter regions 
acts as a repressor; (d) The formation of a structurally different XY-complex (Z2) 
competitively inhibits the DNA binding of X an d Y molecules; fe) Comb ination of 



(c) and (d) as suggested for GATA-1 and PU.l (|Zhamr et al!ll99 9. 2000) 



2.2 Model equations 



With these assumptions one can write down a set of chemical reaction equa- 
tions which underly the system dynamics. 



The processes of specific and unspecific transcription activation (see Fig. 



1(a), (b)) are described by equations (l)-(4). 



x + x + d x ^ d xx ; /;;'• s a /;;:'• • .v (i) 

Y + Y + D X U Df- Df ^ Df + X (2) 

Y + Y + D y 5 Djf; %D™ + Y (3) 

x + x + D„§^ ; "- d xx • y (4) 



Herein we made the simplifying assumption that the DNA binding of X and 
Y always occurs as the binding of homodimers. That means, that the sequen- 
tial binding of two monomers, as the second possibility of DNA binding, is 
not consider. The process of dimerization, as well as the DNA binding and 
dissociation, are regarded to be in quasi steady state. 

Here and throughout the paper, the Ki = fcj/fcj (i = 1, 7) denote the equi- 
librium (dissociation) constants of the reactions, with ki and ki representing 
the forward and backward reaction rate constants, respectively. Finally, it is 
assumed that both transcription factor monomers, X and Y, decay with first 
order kinetics at rates k x and whereas dimer-complexes are assumed to be 
stable. 

The different mutual transcription inhibition mechanisms are illustrated in 
Figs. 1(c)- (e). First, we consider the formation the XF-complex Z2 (see Fig. 
1(d)) 

X + Y^Z 2 . (5) 

Under the quasi steady state assumption Z 2 does not contribute to the math- 
ematical description of the system dynamics. 

As shown in Fig. 1(c), there is also the possibility that X and Y form a 
structurally different heterodimer Zi, which is able to bind to the promoter 
regions, acting as a repressor for X and Y transcription, respectively: 

X + Y + D x ^ D x x v, (6) 
X + Y + D y U D xy . (7) 

As with the promoter binding of X and Y, we collapse dimerization, which 
is assumed to be in quasi steady state, and DNA binding into one process, 
neglecting the sequential binding of monomers. 

Under the posted quasi steady state assumptions, equations (l)-(7) lead to 
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the following set of ordinary differential equations: 



dx _ s x Kix 2 + u x K 2 y 2 , . 

dt ~ *°- X + 1 + K x x 2 + K 2 y 2 + K 6 xy W 
dy = , s y K 3 y 2 + u y K 4 x 2 

dt M 1 + K 3 y 2 + K A x 2 + KTyxy 1 ; 



Details of the derivation are given in Appendix A. 



3 Results 



3.1 Symmetric system 



To analytically derive steady state as well as potential bifurcation conditions, 
we restrict ourself in this section to the special case of a completely symmetric 
system, i.e.: k 0x = k 0y = k , s x = s y = s, u x = u y = u, K x = K 3 , K 2 = K 4 , 
and Kq = K 7 . Using these relations, together with x = ^/K[x, y = y/Kiy, 
k u = K 2 /Ki, k r = Kq/Ki, s = y/K^s/ko, u = y/Kiu/k , and r = k t, the 
system (8), (9) can be written in a dimensionless form as 

dx _ sx 2 + uk u y 2 

dr~ X 1 + x 2 + k u y 2 + A; r xy ' 1 ' 

dy _ sy 2 + uk u x. 2 

dr~ y+ l + A; u x 2 + y 2 + A; r xy' 1 ' 



Equations (10) and (11) are a pair of coupled first order differential equations. 
The steady state (x = y = 0) is defined implicitly by 

sx 2 + w£; u y 2 . , 

(12) 



1 + x 2 + A; u y 2 + A; r xy ' 
sy 2 + ukuX 2 



y ~ l + A; u x 2 + y 2 + A; r xy- 1 } 

The domain of these nullclines for x and y is restricted by the choice of 
parameters as outlined in Appendix B. The intersections of the nullclines 
correspond to the fixed points (x*,y*) of the differential equations (10) and 
(11). Fixed points on the diagonal (x*,x*) are traced under the simplifying 
condition x = y. In this case, equations (12) and (13) can be summarized by 

x* = X * 2 ( S + M ^ (14) 
l+x* 2 {l + k u + k r ) ' { V 
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The first (trivial) fixed point of equation (14) is x* = 0. Having eliminated 
this solution, the remaining quadratic equation yields two further non-trivial 
fixed points at 

(s + uk u ) ± yf (s + uk u ) 2 - 4(1 + k u + k r ) 
X2/3= 2(1 + k u + kr) • (15) 

(xgjXg) and (xg,Xg) are real fixed points on the diagonal for 



s> -uk u + 2^/1 + k u + k r . (16) 

Bifurcation points can be characterized by nullclines intersecting with equal 
slopes. The derivatives of equations (10) and (11) are evaluated to determine 
explicit conditions for bifurcation occurrence on the diagonal, considering s as 
the bifurcation parameter 1 . For simplicity the denominators in equations (12) 
and (13) are defined as P x = l+x 2 + /c u y 2 + /c r xy and P y = l+/c M x 2 +y 2 + A; r xy. 
The partial derivative of equation (12) with respect to y leads to 

, _ (2sxx' + 2uk u y)P x - (sx 2 + uk u y 2 )P' x 

X — p2 { ll > 

x 

with x' = <9x/<9y and P' x = dP x /dy = x'(2x + k r y) + 2A; M y + A; r x. Solving for 
x' yields 

, = 2uk u yP x - (sx 2 + uk u y 2 )(2k u y + Ayx) 
X ~ P 2 - 2sxP, + (sx 2 + uk u y 2 ) (2x + k r y) ' 1 j 

For bifurcation points on the diagonal (x = y), where the denominators P x 
and P y simplify to P* = 1 + x* 2 (l + k u + k r ), equation (18) can be rewritten 

as 

x'(P* 2 -2sxP*+x 3 (s+wA; u )(2 + A; r )) = 2uk u xP* -x 3 (s+uk u )(2k u + k r ) . (19) 

Inserting P* in the form P* = x(s + uk u ) derived from equation (14) and 
neglecting the trivial solution the equality now reads 

x.'(uk u - s + x(2 + kr)) = 2uk u - x(2A; u + k r ) . (20) 



To find the bifurcation points on the diagonal one needs to study the two 
distinct cases for x' = 1 and x' = — 1 (see Appendix C). 

Case I (x' = 1): 

1 Parameter s is chosen to account for changes in the transcriptional activity by 
enhancer actions or modifications in chromatin structure. Furthermore, s is the 
critical parameter that gives rise to the different distinct domains for the nullclines 
as outlined in Appendix B. 
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Equation (20) satisfies the condition x' = 1 at 



s + uk u 

Xx ' =1 " 2(1 + k u + k r ) [Zi) 

This only coincides with the fixed points xj^g derived in equation (15) if the 
expression under the radical in equation (15) vanishes, i.e., x^ = Xg. This is 
true for 

s* = -uk u + 2^1 + k u + k r , (22) 

which corresponds to the condition defined in (16). This implies that the 
"birth" of the fixed points (x^x^) and (xg,Xg) coincides with the bifurca- 
tion condition x' = 1. The sequence of nullclines in Fig. 2(a)-(c) illustrates 
this behavior for an unspecific transcription rate u — 1. The nullclines in 
Fig. 2(a) do not intersect for s < s*, i.e., there is no non-trivial fixed point 
along the diagonal. For s = s\ one common fixed point at x^s*) = Xg(s*) = 
(s* + uk u ) /2(1 + k u + k r ) exists, which marks the bifurcation point depicted 
in Fig. 2(b). For s > s\ two distinct fixed points (x^Xg) and (xg,Xg) exist on 
the diagonal, shown in Fig. 2(c). Whereas the upper point at (x^, x£) is stable, 
the lower one at (xg,Xg) is unstable. The nullclines change qualitatively for a 
further increase in the bifurcation parameter s as shown in the sequence Fig. 
2(d),(e). 

In the case of a smaller unspecific transcription rate u — 0.4 the correspond- 
ing bifurcation is illustrated in Fig. 3(d). The two fixed points (x^x^) and 
(xg, Xg) generated at the diagonal are both unstable as depicted in Fig. 3(e), (f). 
The qualitative differences between the scenarios for small and large unspe- 
cific transcription u are more thoroughly investigated in the subsequent para- 
graphs. 

Case 2 (x' = -1): 

When x' = — 1 equation (20) simplifies to 

3uk u — s 

Xx ' = - X = 2k -2 (23) 

Equating x x / = _i = x^g from equation (15) leads to a dependency on the 
parameters s, u, k u and k r . Two further bifurcation points are obtained at: 



uk u (l + 3Ay + 5A: M ) + 2{k u - 1)^1 - k r - 3k u + IF- 
(-1 + k r + 3k u ) 



ik u (l + 3k r + 5k u ) - 2{k u - 1)^/1 -k r -3k u + 4kl 
(-1 + k r . + 3k u ) 



(24) 
(25) 



To guarantee the existence of these bifurcations at the diagonal, s^ 3 > is 
required. The case S2/3 < s l indicates that bifurcations occur off the diagonal. 
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Figure 2. Deformation of the nullclines for increasing values of the bifurcation pa- 
rameter s. The occurrence of the first bifurcation at s\ is depicted in (b), the second 
bifurcation at s\ = s| in (g). Parameters are k u = 1, k r = and u = 1. The bi- 
furcation parameter is set to s = 1.7 (a), s = = — 1 + 2\/2 pa 1.83 (b), s = 1.9 
(c), s = 1.99 (d), s = 2.01 (e), s = 2.6 (f), s = s* 2 = s* 3 = 3u = 3 (g), s = 3.8 (h). 
Fixed points are marked as follows: trivial fixed point (0,0) - ■. stable/unstable 
fixed point (x^x^) - A/ A, unstable fixed point (x^Xg) - V, stable/unstable fixed 
points off the diagonal - »/o, bifurcation point - □. 
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Figure 3. Deformation of the nullclines for increasing values of the bifurcation pa- 
rameter s for unspecific transcription rate u = 0.4. Parameters are k u = 1 and 
k r = 0. The bifurcation parameter is set to s = 1.9 (a), s = 2.04 (b), s = 2.1 (c), 
s = s* « 2.43 (d), s = 2.8 (e), s = 3.2 (f). Fixed points are marked according to 
the caption in Fig. 2. 

For the special case k u = 1 the conditions for the occurrence of these bifurca- 
tions simplify to = S3 = 3m (given = S3 > si which is true for u > l/\/2). 
Since this condition is valid for both fixed points, (x^Xg) and (x^Xg), it in- 
dicates that the bifurcations occur at the same bifurcation parameter s = 3u. 
Fig. 2(f)-(h) depicts the bifurcations for both fixed points in the case u — 1, 
k u = 1. After a deformation of the nullclines the intersections in Fig. 2(f) still 
represent the fixed points (xg, X2) and (X3, X3) for s < = S3. In Fig. 2(g) the 
nullclines for s = S2 = S3 intersect with the same local slope at X2 as well as at 
Xg. This marks the bifurcation point for both fixed points, that coincides for 
k u — 1. Fig. 2(h) illustrates the new fixed points off the diagonal, which are 
stable bifurcating from (x^x^) and unstable bifurcating from (x^x^). The 
fixed point (x^x^) itself changes the stability and becomes unstable, (x^x^) 
remains unstable as before. 

For small u the condition for the occurrence of further bifurcations S2/3 > s* 
is violated. Numerical results indicate that two saddle-node bifurcations form 
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Figure 4. Bifurcation diagrams x vs. s with k u = 0.8, k r = 0, u = 1 (a) and = 0.4 
(b). The stability of the steady states is coded as follows: solid line - stable, dashed 
line - unstable. 

fixed points off the diagonal at s < s* as depicted in the sequence of nullclines 
in Figs. 3(b), (c). The saddle- node bifurcation on the diagonal is observed at 
s*. For large s these scenarios show a comparable pattern of two up- regulated 
steady states with one high and one low expressed component and a further 
stable fixed point at (0,0) (compare Figs. 2(h) and 3(f)). 

The bifurcation diagrams in Fig. 4 comprise the above findings for k u = 0.8. 
The x-coordinate for the fixed points is shown depending on the bifurcation 
parameter s. For u — 1 (Fig. 4(a)) the birth of two fixed points through a 
saddle-node bifurcation can be seen at s*, given by equation (22). Condition 
(24) defines the occurrence of the pitchfork bifurcation on the upper branch 
(x^x^) at S2, whereas condition (25) is the equivalent for the lower branch 
(x3,Xj|) at S3. Note that the additional condition Sy 3 > s i is fulfilled. The 
upper branch gives rise to three fixed points, one unstable (arising from the 
existing stable fixed point) at the diagonal at Xg and two new stable fixed 
points branching off this axis. For the lower case all three fixed points are 
unstable for s > S3. The inset in Fig. 4(a) enlarges this bifurcation occurring 
at S3. Fig. 4(b) illustrates the equivalent scenario for u = 0.4. The saddle-node 
bifurcation at s* represents the formation of the unstable fixed points (x^x^) 
and (x^Xg). In addition, two further saddle node bifurcations exist that also 
form stable fixed points. Since Sy 3 < s \ these bifurcations do not occur on 
the diagonal. All branches in Fig. 4 that do not represent fixed points on the 
diagonal were determined numerically. 

Fig. 5 provides an overview of regions of multi-stability in the phase space u 
vs. s. Distinct regions with different numbers of stable steady states are iden- 
tified depending on the combination of the dimensionless parameters. Lines of 
separation are determined by equations (22), (24) and, for the lower branch, 
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Figure 5. Phase space diagram u vs. s with k u = 1, k r = 0. The lines separating the 
distinct regions of multi-stability are determined by equations (22), (24) and, for the 
lower branch, numerical results. In the lower left region only one stable fixed point 
at (0, 0) exists. In the region marked with "two stable fixed points" one additional 
up-regulated stable fixed point exists besides the one at (0, 0). In "three stable fixed 
points" region two additional up-regulated stable fixed points exist. The dashed 
horizontal lines correspond to the sequences of nullclines in Figs. 2 and 3. 



numerical results. The sequence of nullclines given in Fig. 2 is illustrated by 
the dashed line at u — 1, with the dots referring to the subfigures for varying 
s. The dashed line at u — 0.4 gives a similar representation, with its corre- 
spondence in Fig. 3. 



Figs. 2 and 3 both indicate that the basin of attraction for the fixed point at 
the origin (0, 0) is separated from the basins of attraction of the up-regulated 
stable states by a set of unstable fixed points. The sequences of graphs also 
illustrate that these unstable fixed points move towards the fixed point at (0, 0) 
for increasing s, thus continously reducing the size of its basin of attraction. 
However, this size characterizes the stability of the fixed point at the origin 
(0, 0) in response to external perturbations. Unlike the intermediate stable 
steady state, arising from the bifurcation at si depicted in Fig. 4(a), where 
a dynamically increasing s inevitably leads to one of the two up-regulated 
fixed points, the escape from the fixed point (0, 0) needs to be triggered by a 
perturbation that exceeds the size of its basin of attraction. Given the position 
of the unstable fixed point at the diagonal (x 3 , Xg) as a function of s in equation 
(15), an appropriate measure for the size of the basin of attraction is provided. 
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3.2 Asymmetric system 



As indicated bv IZhang etafl (|l999l l2000h the inhibition of PU.l by GATA- 
1 and the converse are based on different mechanisms. The formation of the 
PU.l-GATA-1 complex, which we refer to as a .^-complex, prevents free tran- 
scription factors from binding to their specific DNA binding sites. A com- 
petitive inhibition in this form affects both transcription factors, although 
Zhang et al. (2000) do not explicitly outline the consequences of binding of the 



PU.l-GATA-1 complex to the PU.l binding site. On the other hand, GATA- 
1 prevents the binding of c-Jun to the DNA bound PU.l protein and thus 
disables the transcription initiation of PU.l. This process explicitly targets 
the PU.l binding sites and introduces a functional asymmetry of inhibition 
mechanisms. 



The mathematical counterpart of this asymmetry is a specific binding rate 
Kq > while keeping K-j = (see equations (8), (9)). In terms of the di- 
mensionless formulation in equations (10) and (11) this translates into two 
different rate constants k Tx > and k r = 0. The additional binding mode 
(k rx > 0) can be interpreted as a reduction in the transcriptional activity of 
the X gene conferring a disadvantage relative to Y. 

For any k Tx > there is a symmetry breaking which shifts the previously 
observed bifurcations off the diagonal and destroys the pitchfork bifurcation 
observed in the symmetric case for large u. The two up-regulated stable fixed 
points are not created instantaneously by the transformation of a previous 
stable state at the diagonal, but the initial stable point remains unchanged 
while a further (saddle node) bifurcation forms the second up-regulated stable 
point alongside with one unstable fixed point. This scenario is shown in the 
sequence of nullclines in Fig. 6(a)-(c). The parameter k rx regulates the distance 
between the up-regulated stable points and the extension of their basins of 
attraction. This is visualized in the bifurcation diagrams in Fig. 6(d)-(f) for 
different values of k Tx . 

For small unspecific transcription rates u, where in the symmetric case the ad- 
ditional up-regulated stable states are created off the diagonal, no qualitative 
changes are introduced by the functional asymmetry. 

The introduction of asymmetry is not necessarily based on different interaction 
mechanisms. It is plausible that auto-regulative transcription activation does 
not require identical transcription rates for the genes of interest. This can be 
described by relaxing the symmetry assumption of Section 2.1, which leads to 
gene specific transcription rates s x and s y . This asymmetry in transcriptional 
activity results in a qualitatively similar symmetry breaking as in the case 
of the mechanistic asymmetry: the pitchfork bifurcation occurring for large 
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Figure 6. The deformation of the nullclines for the case k Tx = 0.1 is depicted in 
figures (a) to (c). Parameters are k r = 0, u = 0.8, k u = 1, and s = 1.99 (a), 
s = 3.0 (b), s = 6.8 (c). The trivial fixed point at (0,0) is marked by ■, the 
stable/unstable fixed points off the diagonal by •/o. The qualitative change in the 
bifurcation behavior is shown in figures (d) to (f). The bifurcation diagrams are 
shown for increases in the asymmetry parameter k Tx . Parameters are k Ty = 0, u = 1, 
k u = 0.8 and k Tx = 0.0 (d), k Tx = 0.01 (e), and k Tx = 0.1 (f). Solid lines indicate 
stable, dashed lines unstable fixed points. 

u is replaced by a remaining stable state alongside a saddle-node bifurcation 
forming the second up-regulated stable state (data not shown). The magnitude 
of the difference in the specific transcriptions rates s x and s y regulates the 
distance between the up-regulated stable states in the phase plane. 

In a scenario where asymmetry of interaction mechanisms occurs alongside 
an asymmetry in the specific transcription rates, the effects on the system 
behavior combine, either amplifying or compensating each other. 



3.3 Over-expression scenarios 



Induced over-expression of a certain critical component is a common experi- 
mental method to study interaction dynamics between different transcription 
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facto r s and has also been applied to the GATA-1 / PU.l system ([Nerlov et al 



200(1 iRekhtman et all fl999L IZhang et all l2000h . These experiments provide 



insight in the stability of the system, interaction time scales, and the role of 
co-factors and interaction mechanisms. We have applied an over-expression 
impulse of amplitude a oe and duration d oe to the model system given in equa- 
tions (10), (11). Characteristics of the dynamic response are only valid under 
the outlined steady state assumptions. A qualitative overview of the simu- 
lation results is presented in Fig. 7. Starting from a fully symmetric system 
as studied in Section 3.1 where, for large s, the system is in one of the two 
up-regulated states (characterized by one high and one low expressed tran- 
scription factor) two modes of over-expression are applied: a short impulse 
over-expression of the lower expressed component and a long and steady over- 
expression of the same component. Not surprisingly, the model reacts to the 
over-expression with two distinct scenarios, depending on the intensity of the 
impulse. For a subcritical over-expression the system returns to the previous 
expression level (indicated in Figs. 7(a) and (d)), whereas for a supercriti- 
cal situation the former expression state is reversed (indicated in Figs 7(b), 
(c), (e) and (f)). Translating this picture into the x vs. y phase plane, the 
supercritical over-expression corresponds to a change from one basin of at- 
traction to another, induced by a crossing of the separatrix. Most available 
experimental techniques to artificially induce gene expression lead to a mas- 
sive over-expression that significantly exceeds physiological levels, a scenario 
still underestimated by Figs. 7(c) and (f). A sensitively tuned expression ex- 
periment is more promising to elucidate critical intensities and time scales 
necessary to induce a permanent shift in the genetic expression patterns and 
thus to characterize the stability of the initial states. 



4 Discussion 



The presented model of transcription factor interaction is based on principles 

of coupled feedback regulations, whic h have previously been proposed for the 

descr i ption of general genet i c switches dBecskei et alll200ltlCinauin and Demongeotl 
2005; iFrancois and Hakiml . l2004t iGardner et all l2000t iGlass and Ka uffman. 
1973) and the modeling of p r okaryotic gene regu lation ( McAdams and Arkinl . 



1998; ISantillan and Mackevl . Il998l 1200 lL I2004J) . Here, specific experimental 



knowledge of activation and inhibition mechanisms of two transcription fac- 
tors (GATA-1 and PU.l), which play a key role in the myeloid/erythroid dif- 
ferentiation process of hematopoietic progenitor cells, is incorporated in this 
general framework. 

Our model analysis particularly focuses on the investigation of the steady 
states of transcription factor expression and there dependence on parameter 
changes. In this context, we are able to analyze the experimentally suggested 
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Figure 7. Scenarios for sub- and supercritical over-expression. In the subcritical sce- 
narios (a) and (d) the transcription factor concentrations remain at the same fixed 
point, whereas in the supercritical scenarios (b),(c),(e) and (f) the over-expression 
leads to a change of the basin of attraction, resulting in an different final value of 
the transcription factor concentrations. Over-expression is applied as a short term 
influence at time t = 30 ((a),(b) and (c)) and long term influence starting at time 
t = 20 ((d), (e) and (f)). Parameters are s = 5, u = 1, k u = 1, and k r = 0. The 
over-expression is applied with amplitude a oe and duration d oc : a oc = 3 and d oe = 1 
(a), a oc = 4 and d oc = 1 (b), a oc = 8 and d oc = 1 (c), a oe = 0.3 and d oe = 50 (d), 
a oe = 0.32 and d oe = 50 (e), a oe = 2.5 and d oe = 50 (f). 

feedback structures and their effects on the system behavior under various 
conditions. 

To facilitate the mathematical analysis, a number of simplifications have been 
made. We interpret the transcription factors described in the model (X and 
Y) as representatives of a more complex factor formation rather than an ex- 
plicit model of PU.l and GATA-1 alone. Also, we are aware that most of the 
statements resulting from the model analysis are only semi-quantitative in the 
sense that for all model parameters, as there are DNA binding-, decay-, and 
transcription-rates, no experimentally determined estimates are available for 
the investigated system. In the same line of argumentation, details of the tran- 
scription/translation process, like the DNA binding sequence of transcription 
factor molecules and the delay induced by the processes of transcription and 
translation, have been excluded from the analysis. Although such phenomena 
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can in fluence the dynamics of the system ( Bundschuh et al. , 20031 : Vilar et al. 
20021 ). these effects are speculative since detailed information about relevant 



rates and time scales are not available. The simplifications arising from the 
quasi steady state assumption outlined in Section 2 for dimerization and DNA 
binding impose further limitations on o ur model with respect to the exact 
description of the system dynamics (c.f. IPirone and Elstonl . 12004 ). However, 



these simplifications do not effect the steady state behavior, and, thus, do not 
alter the results derived in Section 3. 



The functional role of the so called priming behavior is a question of par- 
ticular biological relevance which is addressed by this model. It has been 
suggested that low level co-expression of multiple transcription factors, spe- 
cific for diff erent lineages , might be a characterist i c of (h e mato poietic) tissue 
stem cells (Akashi, 120051 : ICross and Enver . 1997 : Orkin , 2000j ). However, it 
is currently unclear whether priming corresponds to a stable state of low 
level co-expression or to a truly zero-expression overlaid by some random ex- 
pression noise. Furthermore, there is a hypothesis that lineage specification 
induction might be a two stage process with a primary initialization of tran- 
scription factor network interaction (i.e., a transition from no expression to 
low l evel co-expression) and a s econdary network-induced differentiation pro- 
cess ( Enver and Greavesl . 1998f ). This perspective immediately leads to the 
questions under which conditions such a two stage process can be established 
and whether such a sequence of different activation states of the transcrip- 
tion factor network requires (multiple) external induction signals or whether 
it represents a system inherent development. 

The suggested model generates two characteristic modes of system stability 
depending on the magnitude of the specific transcription rate s: For small s 
only the trivial fixed point (0, 0) exists; for large s two additional up-regulated 
stable states are observed that are marked by the dominance of one factor 
over the other (dominated co-expression). These modes are maintained in- 
dependently of a mechanistic or parametric asymmetry. Assuming a differ- 
entiation initiation b y increasing the transc r iption rate s (e.g. by change s in 
chromatin structure (|Berger and Felsenfeld . 2001 : Rosmarin et all 2005f ) or 
by alterations in activation/inhibition complexes ( Humel . l2000j )). the transi- 
tion between the different stable states is the central mechanism characterizing 
lineage specification. 



Within the proposed biological framework the trivial fixed point at (0,0), 
which exists for all values of s, can be identified with the undifferentiated 
state of a cell where neither activation nor decision processes are observed. It 
should be mentioned that stability of this fixed point is specific for the out- 
lined _jnodeJ_and_lms iwtbeen observed for the general case of a toggle switch 
(c.f. iGardner et all l2000l ). In logical extension, the two up-regulated stable 
fixed points, observed for large s, would be interpreted as expression states 
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promoting one or the other lineage. These distinct states are characterized 
by a high auto-regulative expression of one dominating factor and a reduced 
expression of the antagonistic factor. 

In Sections 3.1 and 3.2 it has been demonstrated that, for increased unspecific 
transcription rates u (but only in this case!), a further stable fixed point exists 
for intermediate s prior to the formation of the two distinct states of domi- 
nated co-expression. This particular fixed point is characterized by a balanced 
low level co-expression of the two antagonistic factors where no final commit- 
ment decision has been made. The resulting transition sequence between three 
distinct regions of multi-stability can be interpreted as a possible explanation 
for a two stage differentiation process mentioned above. 

The induction of a system change from the stable trivial fixed point to the 
dominated or, if existent, balanced low level co-expression state, needs to be 
triggered either by a stochastic background expression or by an active impulse 
on the system. The unstable fixed point separating the trivial from the up- 
regulated stable states is an indicator of the size of the basins of attraction. 
The observation that the unstable fixed point approaches the trivial one for 
increasing s indicates that the magnitude of the perturbation to introduce a 
transition from the zero-state to the co-expression states decreases in the same 
fashion: for a sufficiently large s even a small perturbation is able to initiate 
differentiation. 

Concluding from these results, there are two different scenarios to explain 
the experimentally suggested priming behavior within the proposed model 
framework: (1) Priming might be considered as the existence of perturbations 
in the expression of transcription factors, imposed on a zero-expression state 
represented by the trivial fixed point at (0, 0), either in the form of stochastic 
background fluctuations (functional noise) or by active impulses. In this sce- 
nario, the perturbations are necessary components of the regulatory system 
to induce a differentiation process. It points to the potential role of stochastic 
effects in the context of decision making in stem cell differentiation as fre- 
quently suggested (see iKaern et all (|2005h for a review). (2) In contrast to 



this scenario, priming can also be explained by the balanced low level co- 
expression state, which becomes unstable for increasing specific transcription 
rates. Due to this parameter dependent loss of stability, this scenario would 
lead to differentiation without the need for external perturbations 2 . However, 
the balanced low level co-expression state is only existent if there is a certain 
degree of unspecific transcription. 

Currently, our results do not allow to decide between the two scenarios. The 



2 To be precise: An infinitesimal perturbation is required to escape from the un- 
stable fixed point. Fluctuations of this magnitude are present in any "real world" 
system. 
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introduction of artificial differentiation impulses of different intensities on un- 
committed cells might be an appropriate way to tackle this question exper- 
imentally. Whereas, a low level co-expression priming (like in scenario (2)) 
would be unaffected by these perturbations, the system could be enforced to 
escape the priming status in scenario (1). Moreover the existence and the 
stability of the different stable system states depend sensitively on the model 
parameters. Due to the lack of available data on transcription and binding 
rates, we are currently not able to specify the biological relevant regimes more 
rigorously. Any experimental approximation of binding and transcription rates 
for the involved components supports the identification of the nature of prim- 
ing. 



The over-expression scenarios presented in Section 3.3 fail to explain experi 
mental findings described by several authors ( Nerlov et al. . 2000l ; Rekhtman et al. 



2003t IZhang et all fcoOO). In spite of the induced up-regulation of one tran- 



scription factor it was observed that the transcription level of the antagonis- 
tic transcription factor remained more or less constant. These observations 
are in contrast to the model results presented here, in which the induced 
over-expression of the initially low expressed factor shifts the equilibrium to 
the opposing co-expression state. Retaining our model assumptions, a po- 
tential interpretation can be given as follows: One of the major functions of 
transcriptional regulators like GATA-1 and PU.l is the activation of a set of 
lineage-specific genes which include further transcription and growth factor s 
as well as functional comp onent s of the committe d lineages ([Tenen . 2003). 
In Sieweke and Gral ( 1998j ) and Tsai et al. (1991), the authors point to a 
continuously modulated set of cooperative lineage-inherent transcription fac- 
tors changing with the state of differentiation. Such secondary complexes of 
transcription factors could in turn act as activators of the initial transcription 
factor, substituting for a simple aut o-regulation and thus stabilizing the initial 
up-regulation pattern ( Humel . 2000l ). In such a scenario our model would only 
account for the initial switching process. The experimentally observed stable 
transcription level of the antagonistic factor in over-expression experiments 
could be interpreted as a substitution of the auto-regulation by secondary 
transcription factor complexes. 



Summarizing, the presented model is able to provide a quantitative explana- 
tion for possible mechanisms underlying lineage specification control in eu- 
karyotic systems. It is able to generate parameter dependent changes in the 
system behavior, with alteration of the number of possible stable steady states. 
Specifically, the model explains states of stable co-expression as well as the 
situation characterized by an over-expression of one factor over the other. The 
conditions inducing shifts from one to another stable state (e.g. parameter 
choice, degree of system disturbances), however, depend in a sensitive manner 
on the assumed activation and inhibition mechanisms. Using the mathematical 
model, we were able to test several combinations of experimentally described 
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feedback mechanisms with respect to their influence on the resulting stable 
states and provide possible explanations for the experimentally suggested dif- 
ferentiation priming of stem cells. 
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A Derivation of transcription factor dynamics 



It is assumed that the transcription of transcription factors X and Y requires 
the existence of activator complexes, i.e., the binding of X or Y dimers to 
the promoter regions of X (D x ) and Y (D y ), respectively. As described in 
Section 2.2, we distinguish between a specific (see equations (1),(3)) and an 
unspecific (equations (2), (4)) transcription activation. Furthermore, there is 
the possibility that X and Y can act jointly as a repressor dimer Zi, inhibiting 
the DNA binding of the X and Y activator dimers (see equations (6), (7)). 

The total amount of promoter sites for X and Y can be specified as the sum of 
unbound (free) and occupied (by repressor or activator molecules) promoter 
regions, i.e., 

Dl% = D x/y + + D^ y + . (A.l) 

Using the equilibrium (dissociation) constants 

Y)xx r)yy f)yy Jjxx 

u x v x v y u y 



2 



D x x 2 ' D x y 2) D y y 2) D y x 

T\xy Jjxy 

obtained from assuming equations (l)-(4), (6), (7) to be in a quasi steady state, 
the fraction of promoter sites contributing to active X and Y transcription is 
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given by 



and 



D xx + D yy K lX 2 D x + K 2 y 2 D x 



D™ D x + K x x 2 D x + K 2 yW x + K 6 xyD x 

_ K x x 2 + K 2 y 2 

~ 1 + K x x 2 + K 2 y 2 + K 6 xy 

Df + D x y x _ K 3 y 2 D y + K A x 2 D y 

Df D y + K 3 yW y + K 4 x 2 D y + K 7 xyD y 

K 3 y 2 + K 4 x 2 
1 + K 3 y 2 + K 4 x 2 + K 7 xy 



(A.3) 



(A.4) 



respectively. 

Taking the (first order) decay rates of X and Y into account, one immediately 
obtains equations (8), (9) by writing down the balance equations for X and 
Y. 



B Domain of the nullclines 



Under the equilibrium assumption equation (12) can be solved for: 



k r x. 2 ± Jk 2 x A - Ak u x(u - x)(sx - 1 - x 2 ) 
^ 2(X) = 2k^j 

which describes a set of nullclines. There is an obvious singularity at x = u. 

The solutions yi/ 2 (x) are real for < h(x) = /c^x 4 — 4A; u x(w — x)(sx — 1 — x 2 ) 
with /i(x) defined as the expression under the square root in the previous 
equation. For k r = the roots of h(x) are located at 

= (B.2) 

4 = u (B.4) 
x^ + /f^ (B.5) 

Real roots at X2/4 exist only for s > 2. Fig. B.l shows the function h(x) for 
s = 1.9 and s = 2.1. In the case s < 2 the parameter tt = X3 restricts the 
definition space of the nullclines to x G [0, X3]. For s > 2 three scenarios exist, 
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Figure B.l. The function /i(x) = fc 2 x 4 — 4/c u x(ii — x)(sx— 1 — x 2 ) is shown for k r = 0, 
u = 1, k u = 1 and s = 1.9 (a) and s = 2.1 (b). 

where the singularity at x = u = Xg marks a boundary for distinct intervals 
in the domain: Xg < x^ (with x G [0,Xg] U [x^x^]), x^ < Xg < X4 (with x G 
[0,x£]U[x£,xJ]) shown in Fig. B.l(b), andxj < x^ (with x G [0, x£]U [xj, xj]). 



C Derivation of bifurcation condition 



The nullclines of the symmetric system derived in (12) and (13) are interpreted 
as functions of x and y: 



x = /(x,y) = 

y = #(*,y) 



sx 2 + uk u y 2 



1 + x 2 + /c u y 2 + fc r xy' 

sy 2 + mA: m x 2 
1 + /c u x 2 + y 2 + A; r xy ' 



(C.l) 
(C.2) 



To derive bifurcation conditions one has to determine the point of tangency 
of the nullclines f(x,y), g(x,y) at a steady state (x*,y*), i.e. 



df(x,y) 
dy 



dg{x,y) 
dx 



(C.3) 



(x*,y*) 



Generally, it holds for inverse functions h and k = h' 1 that k'(h(x)) = 
(^'(a;)) -1 . Considering only points at the diagonal x = h(x) = y, it follows 
that h'(x) = (k'(x))^ 1 . Assuming identity of the first order derivatives h' and 
k! at some point x* on the diagonal yields, therefore, (h'(x*)) 2 = (k'(x*)) 2 = 1. 

From these statements, it follows that we have to consider the following equal- 
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ities to find the bifurcation conditions for the symmetric system, restricting 
to symmetric steady states of the form (x*, x*): 



df(x,y) 



dg(x,y) 



1|. 



(C.4) 



dx 



(x*,x*) 



dx 



(x*,x*) 
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